SUBROUTINE G_CONT
!-----------------------------------------------------------------------
!
! Definition von G:  Kontinuitaetsgleichung
!

      use primvar,  only : G, X, XA, MD, MR, MUr, Gnorm
      use physco,   only : z1
      use global,   only : tst, relax
      use config,   only : np, rkind
      use geomvar,  only : S_flux, S_vol, S_volA
      use advecvar, only : rho_adv
      use RBvar,    only : rho_in, rho_out


      implicit none

      integer       :: i
      logical, save :: init = .true.

      if (init) then
         if (.not. relax) then
            write(66,"(a)") "G_cont.f90:   Innenrand: drho/dR = 0 ;  Aussenrand: drho/dR = 0"
         else
            write(66,"(a)") "G_cont.f90:   Zum Relaxieren: Innenrand: rho = rho_in ;  Aussenrand: rho = rho_out"
         end if
         init = .false.
      end if


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

      ! innere Pseudozellen: i=1 & i=2
!~       if ( X(MUr,3) * tst / X(MR,3) > +1.e-30_rkind ) then  ! Einstroemung ueber Innenrand

         G(MD,1)        =  X(MD,1)    - rho_in
         Gnorm(MD,1)    =  X(MD,1)
         G(MD,2)        =  X(MD,2)    - rho_in
         Gnorm(MD,2)    =  X(MD,2)
!~
!~       else
!~
!~          G(MD,1)        =  X(MD,1)    - X(MD,2)
!~          Gnorm(MD,1)    =  X(MD,1)
!~          G(MD,2)        =  X(MD,2)    - X(MD,3)
!~          Gnorm(MD,2)    =  X(MD,2)
!~
!~          rho_in         =  X(MD,3)
!~
!~       end if


      ! aeussere Pseudozellen + ueberzaehlige skalare Zelle: i=np, i=np-1, i=np-2
!~       if ( X(MUr,np-2) * tst / X(MR,np-2) < -1.e-30_rkind ) then  ! Einstroemung ueber Aussenrand

         G(MD,np)       =  X(MD,np)   - rho_out
         Gnorm(MD,np)   =  X(MD,np)
         G(MD,np-1)     =  X(MD,np-1) - rho_out
         Gnorm(MD,np-1) =  X(MD,np-1)
         G(MD,np-2)     =  X(MD,np-2) - rho_out
         Gnorm(MD,np-2) =  X(MD,np-2)

!~       else
!~
!~          G(MD,np)       =  X(MD,np)   - X(MD,np-1)
!~          Gnorm(MD,np)   =  X(MD,np)
!~          G(MD,np-1)     =  X(MD,np-1) - X(MD,np-2)
!~          Gnorm(MD,np-1) =  X(MD,np-1)
!~          G(MD,np-2)     =  X(MD,np-2) - X(MD,np-3)
!~          Gnorm(MD,np-2) =  X(MD,np-2)
!~
!~          rho_out        =  X(MD,np-3)
!~
!~       end if


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

         G(MD,i) = S_vol(i) * X(MD,i) - S_volA(i) * XA(MD,i)                                                         & ! temporal diff.
                 + S_flux(i+1) * rho_adv(i+1) - S_flux(i) * rho_adv(i)                                                 ! advection

         Gnorm(MD,i) = max( &
              abs( S_vol(i) * X(MD,i) ),                                                                             & ! cell content
              abs( S_flux(i+1) * rho_adv(i+1) ), abs( S_flux(i) * rho_adv(i) )                                       & ! advection
                          )

      end do


END SUBROUTINE G_CONT


